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ABSTRACT 

The long-term evolution of supernova remnants (SNRs) in the primordial interstellar medium (ISM) 
with an inhomogeneous structure is calculated to investigate metal enrichment of the primordial gas. 
For this purpose, we have constructed a parallel 3-D hydrodynamics code iircorporating the radiative 
cooling and self-gravity. The self-gravity and radiative cooling develop the inhomogeneous structure of 
the ISM from a small perturbation with a power-law spectrum. The resultant density ranges from 0.5 
cm~'^ to 10^ cm ~^ . Calculations for a supernova (SN) with the progenitor mass of 20 Mq are performed 
as the first step of a series of our study. It is found from the results that a single SN distributes some 
of newly synthesized heavy elements into a dense filament of the ISM with densities ranging from 100 
to lO"' cm^^ depending on where the SN explodes. Thus, the metallicity [Mg/H] of the dense filaments 
polluted by the SN ejecta becomes —2.7 ± 0.5. From these filaments, the first Population II stars will 
form. This value is in accordance with the previous analytical work (Shigeyama & Tsujimoto) with 
accuracy of ~ 0.3 dex. 

Subject headings: galaxies: ISM — hydrodynamics — ISM: abundances — ISM: structure — methods: 
numerical — supernova remnants 



1. INTRODUCTION 



Kee, & Bertschinger (1988). Most of these previous stud- 



Recently, the diversity in the abundance patterns of 
metal-deficient stars in t he Galactic halo has been ob- 
served by several autho rs (McWilliam et al. 1995; Ryanj 



Norria, fc Beers 1D06[). On the other hand, theoretical 



attem pt s- ( |.\udouBO fc Silk 1995 1 ; | S Mge^ 'ama fc Tcujimo ii[ 



1998) have led to the conclusion that star formation in- 
duced by quite few or a single supernova (SN) event results 
in the obs erved diversity in the abundari ce patterns. In 
particular, Shigeyama fc Tsujimoto (1998| ) have suggested 
that star formation occurs in the thin shell produced by a 
SN explosion and that stars of the next generation inherit 
the abundance pattern thereof. This shell was assumed 
to contain all the heavy elements newly synthesized and 
ejected by the SN. As a result, the abundance pattern of 
new stars is determined by combinations of the SN ejecta 
and the interstellar medium (ISM) swept up by the explo- 
sion. In this scenario, the Rayleigh- Taylor instability is 
supposed to happen to destroy the dense heavy element 
layer into numbers of blobs which penetrate into the su- 
pernova remnant (SNR) shell. It is clear that a spherically 
symmetric SNR in a uniform ISM cannot bring its heavy 
element layer into the thin shell formed immediately be- 
hind the shock front. To check whether the abundance 
in the thin shell produced by a SNR reach es the average 
abundance inside the SNR as predicted by Shigeyama & 
Tsujimoto (1998), it is necessary to calculate hydrody- 
namical evolution of the SNR in three dimensions. 

The e volution of a S NR is characterized by several 
phases ( phevalier 1977| ) as the ejecta-dominated (ED) 
phase, Sedov- Taylor (ST) phase, pressure-driven snow- 
plow (PDS) phase, and momentum-conserving snowplow 
(MCS) phase. Each phase has been investigate d by a num- 



ber of authors , e.g., for the first two phases see PYuelove & 
McKee (1999| ) and for the last two phases see |Cioffi, Mc- 



ies have been restricted to spherically symmetric SNRs in 
the uniform ISM. However, the inhomogeneity of an ISM 
affects the evolution of a SNR. A model of SNR evolu- 



tion in an ISM with a density gradient (Hnatyk fc Petruk 



1999) has shown that the shape of the SNR was easily 
deformed to be non-spherical in ~ 10,000 yr after explo- 
sion. However, our objective is to see how a single SN 
event distributes newly synthesized heavy elements in the 
ISM when the SNR stops expanding at t > 1 Myr after 
explosion. The ISM must also evolve on this time-scale 
due to gravity and radiative cooling. Hence, we need to 
follow the evolution of a SNR in the inhomogeneous and 
evolving ISM for more than 1 Myr. 

In this Letter, we investigate the evolution of a SNR 
originated from a Population HI (Pop III) star in the in- 
homogeneous ISM with the primordial abundances by us- 
ing a 3-D hydrodynamics code. Especially, we will focus 
our attention on the correlation of the abundance of heavy 
elements with the gas density in each cell to identify the 
abundance of heavy elements that would be inherited by 
stars of the next generat ion. Consequently, we can test th e 
hypothesis proposed by Shigeyama fc Tsujimoto (1998). 



2. THE METHOD AND INITIAL MODEL 

To model the three-dimensional evolution of a SNR, we 
have constructed a parallel 3-D hydrodynamics code in- 
corporating self-gravity and radiative cooling of the ISM 



with no heav y elements. We adopt the Godunov type (Go- 
dunov 1959| ) scheme to solv e the hydro dynamical equa- 



tions (the HLLC method in [loro 1997| ). The usage of 
the first order scheme reduces not only the computational 
costs of the 3-D calculation by a facto r of two compared 



with higher order schemes like PPM (Colella, & Wood- 



ward 1984) but also the number of boundary values to 
be transferred in parallel version. To obtain a second or- 
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der scheme in time, the Strang-type dimensional spHtting 
( Strang 1968 ) is used. In addition, we have made two 
modif ications to the usual Godunov type scheme. First, to 



cm (this equals to the initial mean density) is ^ 84 %. 
The velocity dispersion is about 3 km s~^. This value will 



propc "ly follow the adiabatic change of thermal energy in 
a pre-shock region, we adopt the m ethod used in a usual 



cosm ological hydrodynamics code (Ryu, Ostriker, Kang 
& Ce: i 1993). Second, the consist ent multi-fluid advec- 
tion method ( Plcwa fc Miiller 1999| ) is implemented since 
we are interested in how heavy elements ejected by a SN 
explosion mix with the ISM. The 3-D Poisson equation is 
solved by the FFT method to calculate the self-gravity of 
the ISM. To solve the energy equation including radiative 
cooling term, wc use an implicit method, i.e., the Newton- 
Raphson method. 

In the early stage of the universe, the gas contains no 
heavy element and the main coolant at temperatures be- 
low 10** K is hydrogen molecules. However, to solve rate 
equations for the formation and destruction of hydrogen 
molecules is not feasible for a simulation like presented 
here. We note that there are som e attempts to take this 
approach (e.g., Abel et al. 1998). Accordingly, through- 
out this Letter, we use the radiative cooling rates for the 
gas with [Fe/H] = — 2 ^ to mimic the cooling function 
of the very early stage of the universe. The cooling rate 
is computed by MAPPINGS III software by R.S. Suther- 
land (MAPP INGS III is the successor of MAPPINGS II 
described in Sutherland fc Dopita 19"93| ). The radiative 
cooling due to hydrogen molecules is not included in MAP- 
PINGS III. In the low temperature range as low as 1000 
K, the adopted cooling rate is slightly lower tha n the cool- 



ing r ate by hydrogen mo lecules (see Figure 2 of Nakasato 
Mori fe Nomoto 2000). We use the equation of states 



incorporating hydrogen and helium in ionization equilib- 
rium. We have assumed that ions and electrons are in 
thermal equilibrium, since the temperature in the SNR is 
almost always below 10^ K during the evolution. 

The ISM is modeled by a periodic square box composed 
of 150"^ zones. The size of the simulation box is set to 80 
pc to cover the maximum size of a SNR. The simulation 
procedure is divided into following two stages. 

The first stage of the simulation is devoted to archiving 
an inhomogeneous ISM structure. Firstly, an initial den- 
sity field (p(r)) for the first stage of the simulation is g ener- 
ated by the COSMICS package ( [Bertschinger 1995| ) with 
a power spectrum of the density expressed as P{k) cx k^^. 
The initial mean number density of hydrogen atoms is set 
to n = 100 cm~'^. We assume that the initial temperature 
is uniform with T = 100 K. Then, we obtain the initial ve- 
locity field (wi(r)) by integrating the equation of motion for 
one dynamical time (tdyn a/I/ (inGp) ^2.7 Myr in the 
present case, where G is the grav itational constan t). We 
use the Zel'dovich approximation dZel'dovich 1970| ) in this 
integration: Ui(r) — ai{r)t^ya, where ai(r) is the accelera- 
tion field generated by p{r) through the Poisson equation. 
Then we have followed the hydrodynamical evolution of 
the box for further 5<dyn- Due to the self-gravity, filamen- 
tary and knotty structures gradually form like in calcula- 
tions of the cosmological structure formation. In Stdyn of 
the evolution, the inhomogeneity develops in the box with 
the density ranging from 0.5 cm~'^ to 10^ cm~'^. The vol- 
ume filling factor of the low density region with rt < 10^ 

^The abundance ratios of the other heavy elements with respect to 



be use d to test the hypothesis of Shigeyama & Tsujimotc 
(199^ ). In the subsequent paper, we will present a de- 



tailed analysis of the 3-D inhomogeneous structure of the 
ISM obtained in the first stage of our simulation. 

In the second stage of our simulation, we model the ST 
and PDS phases of the evolution of a SNR with the same 
code. For the initial condition of the SNR model, we de- 
posit the SN explosion energy {Esn = 10^^ erg) into a cell 
as thermal energy. At the same time, we add the pro- 
genitor mass (20 Mq as a fiducial value) to the same cell 
and treat 1.7 % (0.34 M p, ) of the progenitor mass as Mg 
dTsujimoto ct al. 199^ ). 

In this Letter, we do not intend to simulate the star 
formation processes. Instead, we show the results for two 
cases : the first SNR event occurs in (1) low density region 
or (2) high density region. We consider the region where 
n is lower than 10^ cm"'^ as the low density region and 
the rest of the box as the high density region. In the case 

(1) , we select an arbitrary cell in the low density region 
as a SN site. After the star formation, the star leaves the 
star forming site and moves around under the gravitational 
field of the ISM. Thus, most of SNR events are expected 
to occur in the low density regions due to the large volume 
filling factor. Also in the case (2), we select an arbitrary 
cell in the high density region as a SN site. In the case 

(2) , we treat 20 Mq of the gas in the cell as the progenitor 
mass. 

3. RESULTS 
3.1. case (1): low density regions 

The first SN occurs in a smooth region where the density 
gradient of the ISM is small. This region corresponds to 
a "hole" (or a tenuous region) of the inhomogeneous ISM. 
The number density (n) of the chosen SN site at the end 
of the first stage is ^ 13 cm""^. After adding 20 Mq ejecta 
to the cell, n becomes ~ 5350 cm~'^. The temperature (T) 
of the SN site becomes ^ 1.2 x 10^ K at the beginning of 
the second stage. Then a strongly shocked region is pro- 
duced and the density of the cell decreases very quickly. 
Although the ejecta soon become nearly spherical in < 0.1 
Myr, they are deformed after the collision with the dense 
structure of the ISM. 

To see the evolution of the shape of the ejecta quantita- 
tively, let V the volume of the metal-enriched regions and 
S the surface area of the metal-enriched regions. We define 
the volume filling factor (V/) as V/L^ and the deforma- 
tion factor {Df) as 5/(47r((3V)/(47r))2/3), where L is the 
size of the simulation box. Figure |l| shows the evolutions 
of Vf (upper panel) and Df (lower panel) with the solid 
line for the case (1). In the PDS phase of a spherical SNR 
in a uniform ISM, the shocked volume evolve as oc t^/"^ 



(Cioffi, McKee, & Bertschinger 1988). If we assume the 
Vf evolves as oc f", the exponent a during the first 1 Myr 
is smaller than unity and then increases to almost unity at 
t ~ 3 Myr. It is likely that the ejecta expand faster than 
the shock front in later phases. The fact that Df shown 
in the lower panel is significantly greater than unity indi- 
cates that the shape of the surface of the ejecta becomes 
very deformed from a sphere due to the Rayleigh- Taylor 

Fe have the "primordial" values in Sutherland & Dopita (1993| ) 
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instability. 

Figure ^ shows the mass of the heavy elements at t ~ 3 
Myr as a function of number density (log n) and metallic- 
ity ([Mg/H]). We note that this mass distribution function 
does not change its shape for last 1 Myr. There are two 
peaks in this diagram. One is located at (logn, [Mg/H]) ~ 
(0, —2) and not so pronounced as the other. The total 
amount of the heavy elements involved in this peak is not 
so large. The low density indicates that these heavy ele- 
ments have been staying near the SN site. The other peak 
is located at (logn, [Mg/H]) ~ (2.1, -2.7). A majority of 
the heavy element is involved in this peak. The swept up 



mass estimated from equation (2) of Shigeyama 
119981) becomes (1021)-° °'' x (S/IO)"'/^ 



Tsuji- 



moto 



3.5 



times greater than the original value in this equation. 
Here, the first factor in the left hand side comes from the 
density and the second from the sound speed (or the ve- 
locity dispersion). Thus the [Mg/H] estimated from this 
equation would be ~ —2.9. The analytical work and 3-D 
calculation h ave shown a fairly good agreem ent in this re- 
spect. What Shigeyama & Tsujimoto (199S) did not men- 
tion is that this peak has a finite width. This metallicity 
distribution can be considered to be a consequence of the 
inhomogeneity of the ISM: A part of the ejecta that pene- 
trates into a dense part of the primordial ISM, mixes with 
it, and tends to get a small metallicity. On the other hand, 
a part of the ejecta that contacts with a tenuous part of 
the ISM, mixes with a small amount of ISM thereof, and 
thus gets a larger metallicity. A part of the ejecta that 
never mixes with the ISM will keep its initial metallicity. 

The resultant [Mg/H] ranges -3.2 < [Mg/H] < -2.2 
with a peak at [Mg/H] ~ —2.7. With the resolution of our 
calculations, i.e. dx ^ 0.53 pc, a cell with logn > 2.43 
includes more than one solar mass. Thus the metallicity 
inherited by the next generation stars may be overesti- 
mated unless the mass of these stars is much smaller than 
one solar mass. 

3.2. case (2): high density regions 

In this case, n and T of the SN site at the beginning 
of the second stage are 11,500 cm" and 5.6 x 10^ K, 
respectively. The evolution of Vf and Df is presented in 
Figure |l| with the dotted lines. We can see that Vf is al- 
ways smaller than in case (1) and -D/ is always higher than 
in case (1). This is because the SN occurs in a dense fila- 
ment that produces a large density gradient. Namely, the 
SNR easily expands perpendicular to the direction of the 
filament but is much slowed down along the direction of 
the filament. The surface of the metal-enriched region at 
t ~ 3 Myr is shown in Figure ^ together with the filamen- 
tary structure of the ISM. The SN site is located at the 
center of the simulation box. The yellow surface represents 
the [Mg/H] = —2.6 iso-surface and the orange filamentary 
shape shows the 3-D density structure. The ejecta expand 
into tenuous regions of the ISM and are deformed by the 
filamentary ISM structure. 

As in case (1), the frequency distribution function of 
the mass of Mg at < = 3 Myr with respect to logn and 
[Mg/H] is shown in the right panel of Figure 0. Although 



there are two peaks at (logn, [Mg/H]) ~ (2.3, —2.6) and 
(logn, [Mg/H]) - (4.0, -2.5) also in this case, both of 
the locations and shap es are different from those for case 
(1). Equation (2) of [Shigeyama fc Tsujimoto (1998D 
gives [Mg/H] ^ —2.8 that agrees with this 3-D calcula- 
tion within a factor of 2. 



4. SUMMARY AND DISCUSSION 

In this Letter, we have presented the first results of 3-D 
hydrodynamical calculations of a SNR originated from a 
Pop HI star in the inhomogeneous ISM with the primor- 
dial abundances. We have constructed an efficient parallel 
numerical code to solve the 3-D hydrodynamical equations 
and the Poisson equations. Using our code, we first evolve 
a random density field and obtain the inhomogeneous ISM 
model. Then, we follow the long-term evolution of a SNR 
to see how a single SN distributes the newly synthesized 
heavy elements into the ISM. 

For the site of the first SN event in the simulation, 
we consider low density and high density environments. 
In both strong shock is produced and then the 

shock front is slowed down by collisions with the filamen- 
tary structure of the ISM. Our calculations indicate that 
a majority of the newly synthesized heavy elements is dis- 
tributed into the gas in filamentary structures to get the 
metaUicity of -3.2 < [Mg/H] < -2.2 depending on the 
SN site. Thus, Pop II stars born from a high density fil- 
ament polluted with the ejecta of a SN will have [Mg/H] 
^ —2.7 on their surfaces. We have compared the values of 
[Mg/H] in dense regions of the simulat ion box with those 



estimated from the analytical work by Shigeyama fc Tsu-^ 
[jimoto (1998 ) and found that they are in agreement with 
each other with an accuracy to 0.3 dex, though the as- 
sumcd density structure of the ISM in Shigeyama & Tsu- 



jimoto (1998| ) is quite different from that in the present 
work. 

To understand metal enrichment process at the early 
stages of galaxy evolution, wc will need to eliminate (or 
identify, at least) the infiuences from numerical artifacts 
that may be caused by the limited resolution of our model. 
First, we needed to deposit thermal energy into one cell 
to initiate a SN. A significant fraction of the explosion en- 
ergy of a real SN with the size of 0.5 pc must be in the 
form of kinetic energy. Thus the ejecta of a real SN are ex- 
pected to penetrate deeper into the filamentary structure 
to get lower metallicity than shown in this Letter. Second, 
the inhomogeneous ISM model obtained in the first stage 
of our simulation procedure may be changed depending 
on a number of model assumptions which include the ini- 
tial power spectrum of the model, the cooling rates due to 
molecular hydrogen, neglect of the dark matter potential, 
the effect of radiative transfer, and etc.. 

We will report results of a detailed study of the inhomo- 
geneous ISM and longer ev olution of the metallicity distr i- 
bution in a separate paper (Nakasato & Shigeyama 2000). 
Calculations for SNe with different progenitor masses will 
be also performed to see the abundance patterns of Pop II 
stars. 
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Fig. 1. — Upper panel : The evolution of the volume filling factor {Vj) of the metal-enriched region. Lower panel : The evolution of the 
deformation factor (Dj) of the metal-enriched region. In both panel, the solid and dotted lines correspond to the case (1) and (2), respectively. 
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Fig. 2. — The mass of the heavy elements at t = 3 Myr with number density (logn) and given metallicity ([Mg/H]). The left and right 
panels correspond to case (1) and case (2), respectively. 
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Fig. 3. — The metal-enriched region at t = 3 Myr for case (2) is rendered by a yellow shaded iso-surface. The surface value corresponds to 
[Mg/H] ~ —2.6. The orange filamentary shape is a volume rendering of the dense structure at the same time. The SN site is located at the 
center of the simulation box shown with the white lines. 



